Spatial data are used across a wide range of fields to support decision-making, including environment, public health, ecology, agriculture, urban planning, economy, and society.
1.1 Areal data
Areal data usually arise when the number of events corresponding to some variable of interest are aggregated in areas.
library(sf)
Warning: package 'sf' was built under R version 4.4.3
Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
Example of areal data. Household income in $1000 USD in neighborhoods in Columbus, Ohio, in 1980.
library(terra)
Warning: package 'terra' was built under R version 4.4.3
terra 1.8.60
d <-rast(system.file("ex/elev.tif", package ="terra"))plot(d)
Example of areal data. Elevation at raster grid cells covering Luxembourg.
1.2 Geostatistical data
For example, we can use air pollution measurements at a number of monitoring stations to predict air pollution at other locations taking into account spatial autocorrelation and other factors that are known to predict the outcome of interest.
Example of geostatistical data. Price per square meter of a set of apartments in Athens, Greece, in 2017.
library(malariaAtlas)
Warning: package 'malariaAtlas' was built under R version 4.4.3
Registered S3 method overwritten by 'malariaAtlas':
method from
autoplot.default ggplot2
d <-getPR(country ="Zimbabwe", species ="BOTH")
Loading ISO 19139 XML schemas...
Loading ISO 19115-3 XML schemas...
Loading ISO 19139 codelists...
Please Note: Because you did not provide a version, by default the version being used is 202406 (This is the most recent version of PR data. To see other version options use function listPRPointVersions)
Warning: The `sourcedata` argument of `isAvailable_pr()` is deprecated as of
malariaAtlas 1.6.0.
ℹ The argument 'sourcedata' has been deprecated. It will be removed in the next
version. It has no meaning.
ℹ The deprecated feature was likely used in the malariaAtlas package.
Please report the issue at
<https://github.com/malaria-atlas-project/malariaAtlas/issues>.
Please Note: Because you did not provide a version, by default the version being used is 202406 (This is the most recent version of PR data. To see other version options use function listPRPointVersions)
Confirming availability of PR data for: Zimbabwe...
PR points are available for Zimbabwe.
Attempting to download PR point data for Zimbabwe ...
Data downloaded for Zimbabwe.
NOTE: All available data for this query was downloaded for both species,
but there are no PR points for P. vivax in this region in the MAP database.
To check endemicity patterns or to contribute data, visit malariaatlas.org OR email us at map@bdi.ox.ac.uk.
ggplot2::autoplot(d)
Warning: The `format` argument of `getShp()` is deprecated as of malariaAtlas 1.6.0.
ℹ The argument 'format' has been deprecated. It will be removed in the next
version. Admin boundaries will be correctly plotted using autoplot without
the argument.
ℹ The deprecated feature was likely used in the malariaAtlas package.
Please report the issue at
<https://github.com/malaria-atlas-project/malariaAtlas/issues>.
Please Note: Because you did not provide a version, by default the version being used is 202403 (This is the most recent version of admin unit shape data. To see other version options use function listShpVersions)
Example of geostatistical data. Malaria prevalence at specific locations in Zimbabwe.
1.3 Point patterns
Point patterns arise when the variable to be analyzed corresponds to the location of events.
library(spatstat)
Warning: package 'spatstat' was built under R version 4.4.2
Loading required package: spatstat.data
Warning: package 'spatstat.data' was built under R version 4.4.2
Loading required package: spatstat.univar
Warning: package 'spatstat.univar' was built under R version 4.4.2
spatstat.univar 3.1-1
Loading required package: spatstat.geom
Warning: package 'spatstat.geom' was built under R version 4.4.2
spatstat.geom 3.3-5
Attaching package: 'spatstat.geom'
The following objects are masked from 'package:terra':
area, delaunay, is.empty, rescale, rotate, shift, where.max,
where.min
Loading required package: spatstat.random
Warning: package 'spatstat.random' was built under R version 4.4.2
spatstat.random 3.3-2
Loading required package: spatstat.explore
Warning: package 'spatstat.explore' was built under R version 4.4.2
Loading required package: nlme
spatstat.explore 3.3-4
Loading required package: spatstat.model
Warning: package 'spatstat.model' was built under R version 4.4.2
Loading required package: rpart
spatstat.model 3.3-4
Loading required package: spatstat.linnet
Warning: package 'spatstat.linnet' was built under R version 4.4.2
spatstat.linnet 3.2-5
spatstat 3.3-1
For an introduction to spatstat, type 'beginner'
plot(clmfires, use.marks =FALSE, pch =".")
Data clmfires is a marked point pattern containing information of each fire, depicts the location of the fires without the mark.
library(spatstat)plot(hamster)
This figure also shows the positions of cell nuclei in a histological section of a tissue from a lymphoma in the kidney of a hamster from spatstat. The nuclei are classified as either “pyknotic” (corresponding to dying cells) or “dividing” (corresponding to cells arrested in the act of dividing).
library(sparr)
Warning: package 'sparr' was built under R version 4.4.3
Welcome to
_____ ___ ____ ____ ____
/ ___// _ \/ _ \/ __ \/ __ \
\__ \/ ___/ __ / ___/ ___/
___/ / / / / / / /\ \/ /\ \
/____/_/ /_/ /_/_/ \__/ \_\ v2.3-16
- type news(package="sparr") for an overview
- type help("sparr") for documentation
- type citation("sparr") for how to cite
data(pbc)plot(unmark(pbc[which(pbc$marks =="case"), ]), main ="cases")axis(1)axis(2)title(xlab ="Easting", ylab ="Northing")
plot(unmark(pbc[which(pbc$marks =="control"), ]), pch =3, main ="controls")axis(1)axis(2)title(xlab ="Easting", ylab ="Northing")
Examples of point patterns. Top: Locations of fires in Castilla-La Mancha, Spain, between 1998 and 2007. Bottom: Locations and types of cells in a tissue.
1.4 Spatio-temporal data
Spatio-temporal data arise when information is both spatially and temporally referenced. Thus, we can consider spatial data as temporal aggregations or temporal snapshots of a spatio-temporal process.
1.5 Spatial functional data
Spatial functional data arise when the three types of spatial data (areal, geostatistical, and point patterns) are combined with random functions.
1.6 Mobility data
Besides the three classical types of spatial data (i.e., areal, geostatistical, and point patterns), we can also consider other spatial data such as flows containing the number of individuals or other elements moving between locations.
library("epiflows")
Warning: package 'epiflows' was built under R version 4.4.3
epiflows is loaded with the following global variables in `global_vars()`:
coordinates, pop_size, duration_stay, first_date, last_date, num_cases
Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` (with slightly different semantics) to convert to a
tibble, or `as.data.frame()` to convert to a data frame.
ℹ The deprecated feature was likely used in the epiflows package.
Please report the issue at <https://github.com/reconhub/epiflows/issues>.
# A tibble: 15 × 8
id location_population num_cases_time_window first_date_cases
<chr> <dbl> <dbl> <fct>
1 Espirito Santo 3973697 2600 2017-01-04
2 Minas Gerais 20997560 4870 2016-12-19
3 Rio de Janeiro 16635996 170 2017-02-19
4 Sao Paulo 44749699 200 2016-12-17
5 Southeast Brazil 86356952 7840 2016-12-17
6 Argentina NA NA <NA>
7 Chile NA NA <NA>
8 Germany NA NA <NA>
9 Italy NA NA <NA>
10 Paraguay NA NA <NA>
11 Portugal NA NA <NA>
12 Spain NA NA <NA>
13 United Kingdom NA NA <NA>
14 United States of … NA NA <NA>
15 Uruguay NA NA <NA>
# ℹ 4 more variables: last_date_cases <fct>, length_of_stay <dbl>, lon <dbl>,
# lat <dbl>
// flows
# A tibble: 100 × 3
from to n
<chr> <chr> <dbl>
1 Espirito Santo Italy 2828.
2 Minas Gerais Italy 15714.
3 Rio de Janeiro Italy 8164.
4 Sao Paulo Italy 34039.
5 Southeast Brazil Italy 76282.
6 Espirito Santo Spain 3270.
7 Minas Gerais Spain 18176.
8 Rio de Janeiro Spain 9443.
9 Sao Paulo Spain 39371.
10 Southeast Brazil Spain 88231.
# ℹ 90 more rows
We can use this data to create an epiflows object called ef that allows us to use the prediction and visualization functions.
2 Spatial data in R
2.1 Vector data
The sf package allows us to work with vector data which is used to represent points, lines, and polygons.
2.1.1 Shapefile
Vector data are commonly stored in the shapefile format.
A shapefile is not a single file, but a collection of related files that work together.
The three mandatory components are:
.shp – contains the geometry data,
.shx – a positional index that allows quick access to the geometry,
.dbf – stores attribute data for each feature.
Additional optional files may include:
.prj – specifies the map projection,
.sbn and .sbx – spatial index files,
.shp.xml – contains spatial metadata in XML format.
Map of the first attribute of the ’sf‘ object representing the counties of North Carolina, USA.
2.2 Raster data
Raster data (also referred to as grid data) is a spatial data structure that divides the region of study into rectangles of the same size called cells or pixels, and that can store one or more values for each of these cells.
2.2.1 GeoTIFF
Raster data often come in GeoTIFF format which has extension .tif.
class : SpatRaster
size : 90, 95, 1 (nrow, ncol, nlyr)
resolution : 0.008333333, 0.008333333 (x, y)
extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : elev.tif
name : elevation
min value : 141
max value : 547
plot(r)
Map of raster data representing the elevation of Luxembourg.
2.3 Coordinate Reference System
2.3.1 Geographic CRS
In a geographic CRS, latitude and longitude values are used to identify locations on the Earth’s three-dimensional ellipsoid surface.
Latitude values range from –90° at the South Pole to +90° at the North Pole. Longitude values measure the angle east or west of the Prime Meridian.
Latitude and longitude coordinates may be expressed in degrees, minutes, and seconds, or in decimal degrees.
Projected CRS
Projected CRSs use Cartesian coordinates to reference a location on a two-dimensional representation of the Earth.
World maps with Mercator (left) and Robinson (right) projections.
A common projection is the Universal Transverse Mercator (UTM) projection. This projection is conformal, meaning it preserves angles and therefore maintains shapes over small areas.
2.3.3 EPSG codes
Most common CRSs can be specified using their EPSG (European Petroleum Survey Group) codes or their Proj4 strings. Common spatial projections can be found at https://spatialreference.org/ref/. Details of a given projection can be inspected using the st_crs() function of the sf package. For example, EPSG code 4326 refers to the WGS84 longitude/latitude projection.
st_crs("EPSG:4326")$Name
[1] "WGS 84"
st_crs("EPSG:4326")$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"
st_crs("EPSG:4326")$epsg
[1] 4326
2.3.4 Transforming CRS with sf and terra
library(sf)pathshp <-system.file("shape/nc.shp", package ="sf")map <-st_read(pathshp, quiet =TRUE)# Get CRS# st_crs(map)# Transform CRSmap2 <-st_transform(map, crs ="EPSG:4326")# Get CRS# st_crs(map2)
library(terra)pathraster <-system.file("ex/elev.tif", package ="terra")r <-rast(pathraster)# Get CRS# crs(r)# Transform CRSr2 <- terra::project(r, "EPSG:2169")# Get CRS# crs(r2)
Alternatively, if we want the transformed data to align exactly with other raster data we are using, we can project it using an existing raster that has the desired geometry.
# x is existing raster# r is raster we project# r2 <- terra::project(r, x)
2.4 Old spatial packages
Before the sf package was developed, the sp package was used to represent and work with vector spatial data. The sp, rgdal (Bivand et al., 2023), rgeos (Bivand and Rundel, 2022), and maptools (Bivand and Lewin-Koh, 2022) packages are no longer maintained and have been retired.
A sf object contains the following objects of class sf, sfc and sfg: • sf (simple feature): each row of the data.frame is a single simple feature consisting of attributes and geometry. • sfc (simple feature geometry list-column): the geometry column of the data.frame is a list-column of class sfc with the geometry of each simple feature. • sfg (simple feature geometry): each of the rows of the sfc list-column corresponds to the simple feature geometry (sfg) of a single simple feature.
Simple feature geometries (sfg objects) can be, for example, of type POINT (single point), MULTIPOINT (set of points), or POLYGON (polygon), and can be created using the functions st_point(), st_multipoint(), and st_polygon(), respectively. Here, we create an sf object containing two single points, a set of points, and a polygon, each with one attribute.
library(sf)library(ggplot2)# Single points (point as a vector)p1_sfg <-st_point(c(2, 2))p2_sfg <-st_point(c(2.5, 3))# Set of points (points as a matrix)p <-rbind(c(6, 2),c(6.1, 2.6),c(6.8, 2.5),c(6.2, 1.5),c(6.8, 1.8))mp_sfg <-st_multipoint(p)# Polygon: sequence of points forming a closed, non-self-intersecting ring# The first ring denotes the exterior ring;# zero or more subsequent rings denote holes within itp1 <-rbind(c(10, 0), c(11, 0), c(13, 2),c(12, 4), c(11, 4), c(10, 0))p2 <-rbind(c(11, 1), c(11, 2), c(12, 2), c(11, 1))pol_sfg <-st_polygon(list(p1, p2))# Create sf objectp_sfc <-st_sfc(p1_sfg, p2_sfg, mp_sfg, pol_sfg)df <-data.frame(v1 =c("A", "B", "C", "D"))p_sf <-st_sf(df, geometry = p_sfc)# Plot single points, multipoint, and polygonggplot(p_sf) +geom_sf(aes(col = v1), size =3) +theme_bw()
sf object representing two single points, a set of points, and a polygon, with one attribute.
3.3 st_*() functions
Common functions to manipulate sf objects include:
• st_read() – reads an sf object
• st_write() – writes an sf object
• st_crs() – gets or sets a new coordinate reference system (CRS)
• st_transform() – transforms data to a new CRS
• st_intersection() – intersects sf objects
• st_union() – combines several sf objects into one
• st_simplify() – simplifies an sf object
• st_coordinates() – retrieves coordinates of an sf object
• st_as_sf() – converts another object type to an sf object
although coordinates are longitude/latitude, st_union assumes that they are
planar
The boundaries of a map can be simplified with the st_simplify() function.
# Simplifymap_utm <-st_transform(map, 32617)map_simplified <-st_simplify(map_utm, dTolerance =10000)ggplot(map_simplified) +geom_sf(fill ="lightblue", color ="gray40") +theme_minimal() +labs(title ="Simplified North Carolina Map (10 km tolerance)")
3.4 Transforming point data to a sf object
The st_as_sf() function allows us to convert a foreign object into an sf object. For example, we can have a data frame containing the coordinates of several locations along with their attributes, which we can then convert into an sf object.
We can use the st_intersects() function from sf to count the number of points within the polygons of an sf object. The returned object is a list containing the feature IDs that intersect each polygon. We can then use the lengths() function to calculate the number of points inside each feature.
library(sf)library(ggplot2)# Map with divisions (sf object)map <-read_sf(system.file("shape/nc.shp", package ="sf"))# Generate random points within the polygons (sf geometry list-column sfc)points <-st_sample(map, size =100, exact =FALSE)
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
# Convert points (sfc) into an sf object so it can be used in ggplot easilypoints_sf <-st_sf(geometry = points)# Plot map and pointsggplot()+geom_sf(data= map)+geom_sf(data=points)
Points within polygons.
# Intersection: which points fall in which polygonsinter <-st_intersects(map, points_sf)
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
# Add point count to each polygonmap$count <-lengths(inter)# Plot number of points within each polygonggplot(map)+geom_sf(aes(fill = count))
Number of points within polygons.
3.6 Identifying polygons containing points
library(sf)library(ggplot2)library(dplyr)
Attaching package: 'dplyr'
The following object is masked from 'package:nlme':
collapse
The following objects are masked from 'package:terra':
intersect, union
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
# Map with divisions (sf object)map <-read_sf(system.file("shape/nc.shp", package ="sf"))# Generate random points within the polygonspoints <-st_sample(map, size =3) %>%st_as_sf()
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
# Find which polygon each point falls inside# (note: first argument = points, second = polygons)inter <-st_intersects(points, map)
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
# Add column with the name of the area containing each pointpoints$areaname <- map[unlist(inter), "NAME", drop =TRUE]# Show points with corresponding area namepoints
Simple feature collection with 3 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -79.0806 ymin: 34.94656 xmax: -76.03757 ymax: 35.96621
Geodetic CRS: NAD27
x areaname
1 POINT (-79.0806 34.94656) Hoke
2 POINT (-76.03757 35.91925) Tyrrell
3 POINT (-78.32428 35.96621) Franklin
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data
3.7 Joining map and data
Sometimes, a map and its corresponding data are available separately, and we may want to create an sf object that combines the map with its related data for easier manipulation and plotting. This can be done by joining the map and data using the left_join() function from the dplyr package.
library(rnaturalearth)
Warning: package 'rnaturalearth' was built under R version 4.4.3
map <-ne_countries(returnclass ="sf")
library(wbstats)
Warning: package 'wbstats' was built under R version 4.4.3
Then, we use the left_join() function from dplyr to merge the map and data, specifying the variable(s) to join by using the by argument. In this example, we use the ISO3 standard country codes instead of country names, as names may differ between the map and the data frame. Figure 3.7 shows the resulting map created using ggplot2.
PM2.5 values in each of the world countries in 2016.
Note that when using left_join(), the resulting object’s class matches that of the first argument.
Thus, left_join(sf_object, data.frame_object) produces an sf object, while left_join(data.frame_object, sf_object) produces a data frame.
map1 <-left_join(map, d, by =c("iso_a3"="iso3c"))class(map1)
[1] "sf" "data.frame"
d1 <-left_join(d, map, by =c("iso3c"="iso_a3"))class(d1)
We use the terra functions centroids() to get the centroids of the division polygons, and crds() to obtain their coordinates.
points <-crds(centroids(v))
plot(r)plot(v, add =TRUE)points(points)
Elevation raster, and division and centroids of polygons in Luxembourg.
We can obtain raster values at points using extract(), with the first argument as the SpatRaster object and the second argument as the data frame containing the points.
# data frame with the coordinatespoints <-as.data.frame(points)valuesatpoints <-extract(r, points)cbind(points, valuesatpoints)
# Extracted raster cells and percentage of area covered within each polygonhead(extract(r, v, na.rm =TRUE, weights =TRUE))
ID elevation weight
1 1 NA 0.04545454
2 1 NA 0.10909091
3 1 529 0.24545454
4 1 542 0.46363635
5 1 547 0.68181816
6 1 535 0.11818181
# Average raster values by polygonv$avg <-extract(r, v, mean, na.rm =TRUE)$elevation# Area-weighted average raster values by polygon (weights = TRUE)v$weightedavg <-extract(r, v, mean, na.rm =TRUE, weights =TRUE)$elevation
library(ggplot2)library(tidyterra)
Warning: package 'tidyterra' was built under R version 4.4.3
Registered S3 method overwritten by 'tidyterra':
method from
autoplot.SpatRaster malariaAtlas
Attaching package: 'tidyterra'
The following object is masked from 'package:stats':
filter
# Plot average raster values within polygonsggplot(data = v) +geom_spatvector(aes(fill = avg)) +scale_fill_terrain_c()
# Plot area-weighted average raster values within polygonsggplot(data = v) +geom_spatvector(aes(fill = weightedavg)) +scale_fill_terrain_c()
Average and area-weighted average of elevation values in each of the divisions of Luxembourg.
5 Making maps with R
Maps allow us to effectively communicate spatial information. Here, we demonstrate how to create both static and interactive maps using several mapping packages, including ggplot2 (Wickham et al., 2022a), leaflet (Cheng et al., 2022a), mapview (Appelhans et al., 2022), and tmap (Tennekes, 2022).
The following object is masked from 'package:ggplot2':
last_plot
The following object is masked from 'package:stats':
filter
The following object is masked from 'package:graphics':
layout
g <-ggplot(d) +geom_sf(aes(fill = vble))ggplotly(g)
The sf object passed to leaflet() must have a geographic coordinate reference system (CRS) using latitude and longitude (EPSG:4326). Here, we use the st_transform() function from sf to convert the data d, which has CRS EPSG:4267, to CRS EPSG:4326.
st_crs(d)$epsg
[1] 4267
d <-st_transform(d, 4326)
library(leaflet)
Warning: package 'leaflet' was built under R version 4.4.2
We can also use the addMiniMap() function to include an inset map.
l %>%addMiniMap()
For example, to save the leaflet map l as .png we can proceed as follows:
# webshot::install_phantomjs()# Saves map.html# library(htmlwidgets)# saveWidget(widget = l, file = "map.html")# Takes a screenshot of the map.html created and saves it as map.png# library(webshot)# webshot::install_phantomjs()# webshot(url = "map.html", file = "map.png")
5.2 mapview
library(mapview)mapview(d, zcol ="vble")
Maps created with mapview can also be customized to include elements such as legends and background maps.
The tmap package (Tennekes, 2022) allows us to create static and interactive maps composed of multiple shapes and layers, with various styles.
library(tmap)
Warning: package 'tmap' was built under R version 4.4.2
tmap_mode("plot")
ℹ tmap mode set to "plot".
tm_shape(d) +tm_polygons("vble")
5.7 Maps of point data
library(maps)
Warning: package 'maps' was built under R version 4.4.2
Attaching package: 'maps'
The following object is masked from 'package:viridis':
unemp
d <- world.cities# Select South Africad <- d[which(d$country.etc =="South Africa"), ]# Transform data to sf objectd <-st_as_sf(d, coords =c("long", "lat"))# Assign CRSst_crs(d) <-4326
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_dots()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).[v3->v4] `tm_tm_dots()`: migrate the argument(s) related to the scale of the
visual variable `size` namely 'scale' (rename to 'values.scale') to size.scale
= tm_scale_continuous(<HERE>).
ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
a list: 'size.scale = list(<scale1>, <scale2>, ...)'
5.8 Maps of raster data
library(terra)filename <-system.file("ex/elev.tif", package ="terra")r <-rast(filename)# Transform data to sf objectd <-st_as_sf(as.data.frame(r, xy =TRUE), coords =c("x", "y"))# Assign CRSst_crs(d) <-4326# Plotggplot(d) +geom_sf() +geom_raster(data =as.data.frame(r, xy =TRUE),aes(x = x, y = y, fill = elevation))
To use the leaflet and mapview packages, we transform the data from class terra to RasterLayer with the raster::brick() function.
library(raster)
Attaching package: 'raster'
The following object is masked from 'package:plotly':
select
The following object is masked from 'package:tidyterra':
select
The following object is masked from 'package:dplyr':
select
The following object is masked from 'package:nlme':
getData
[v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
visual variable `col` namely 'palette' (rename to 'values') to col.scale =
tm_scale(<HERE>).
[v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
"brewer.yl_or_rd"
Finally, we call the flowmapblue() function, passing the locations and flows data frames, along with the Mapbox access token, and specifying options such as clustering and animation. If mapboxAccessToken is set to NULL, the flows between locations will still be visualized, but without a background map.
Warning: package 'elevatr' was built under R version 4.4.3
elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'. Use
of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being
deprecated; however, get_elev_raster continues to return a RasterLayer. This
will be dropped in future versions, so please plan accordingly.
library(terra)map <-ne_countries(type ="countries", country ="Switzerland", scale ="medium", returnclass ="sf")d <-get_elev_raster(locations = map, z =9, clip ="locations")
For example, we can obtain the hospitals in Barcelona by specifying its bounding box placebb and using add_osm_feature() with key = "amenity" and value = "hospital" as follows.
hospitals <- placebb %>%opq() %>%add_osm_feature(key ="amenity", value ="hospital") %>%osmdata_sf()
Motorways can be retrieved using key = "highway" and value = "motorway".
motorways <- placebb %>%opq() %>%add_osm_feature(key ="highway", value ="motorway") %>%osmdata_sf()